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ABSTBACT 

Since widespread use of the finite elemenr method began 
in the early 1960's, much effort has been devoted to the 
development of the method itself, while only recently has 
there been any research directed ar minimizing the discreti- 
zation error by a proper selection of ~he element grid. This 
paper examines some recently proposed grid optimization 
techniques and applies them to some one- and two-dimensional 
linear self-adjoint boundary value problems. Guidelines 
requiring minimal software modification are recommended to 
assist the analyst in obtaining improved finite element 
solutions. 
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IN TROPaCTI DH 



The critical concern in any finite elemenT: analysis is 

the adequacy of the selected mesh to provide reliable solu- 
tion results within some reasonable cost range. The goal of 
finite element grid optimization nhen becomes one of 
obtaining maximum solution accuracy for a prescribed anal- 
ysis cost. While this objective is generally not realized in 
today's widespread use of finite element analysis, the effi- 
cient computation of solutions with optimal accuracy will 
become paramount to the engineer as finite element methods 
are applied to increasingly difficult dynamic, nonlinear, 
and evolution problems. 

A. HISTORICAL BACKGROOND 

In the early 1960's, with the help of the high speed 
digital computer, finite element methods revolutionized 
problem solving in engineering. Since that time the major 
research efforts have concentrated on expanding the theoret- 
ical basis of the method and extending its application in a 
variety of fields. Only recently has there been significant 
attention directed at minimizing finite element solution 
errors by properly defining the element grid. Early attempts 
at distributing the nodes and choosing the elements to 
ensure seme degree of confidence in the solution accuracy 
were predominantly dependent upon the analyst's engineering 
judgement and experience, since there were no established 
procedures for accomplishing this objective. Sven these 
attempts towards grid optimization have become largely 
ignored with the advent of automatic mesh generators, which 
have drastically reduced preprocessing costs while 
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accomplishing little in improving solution accuracy. These 
programs automatically construct the element grid, usually 
in a uniform manner, after the user merely defines the 
problem and specifies the number of elements to be used. To 
establish convergence and achieve the desired solution accu- 
racy, the user simply repeats the analysis using a finer 
mesh of uniformly distributed elements while monitoring such 
convergence indicators as successive solution values at 
common nodal points or the asymptotic increase in the energy 
content of the mesh. The often disastrous flaw in such a 
practice is that for nearly degenerate problems which 
exhibit very large gradients, the asymptotic range is only 
entered for an extremely large number of degrees-of- freedom, 
often beyond computer limitations [Sef. 1]. In this case, 
uniform mesh refinement may produce apparent convergence, 
when in fact the solution accuracy is poor. Therefore, the 
need for a finite element grid optimization procedure is 
clearly evident. 

The first formal attempts at finire element grid optimi- 
zation did not begin until the early 1970 's. These early 
approaches involved the inclusion of the nodal coordinates 
as dependent variables in the minimization of the potential 
energy functional [Ref. 2], Unfortunately, the resulting 
system of equations is highly nonlinear and the computa- 
tional effort involved in its solution is sc great that 
similar accuracy can be obtained at a fraction of the 
expense, simply by employing a very fine mesh. Clearly, this 
method does not approach the finite element grid optimiza- 
tion goal of achieving a specified solution accuracy for a 
minimum analysis cost. For this reason, virtually all of the 
grid optimization techniques since developed are based on a 
"near-optimum” strategy whereby nearly-optimal solution 
results are obtained without the computational inefficiency 
of a numerical optimization analysis. The growing emphasis 
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has been on adaptivity, a procedure for efficient construc- 
tion of near-optimum grids by the iterative application of 
some criterion, based on data already computed from the 
solution for a previous grid. This procedure is far more 
efficient than the conventional approach of repeating the 
analysis using successively finer uniform grids. 
Experimental self-adaptive finite element codes have 
recently been developed. Starting from the user’s initial 
idealization, these programs automatically generate a near- 
optimum grid and solve the resulting equations. 

B. IHVESTIGATIVE APPROACH 

In undertaking any numerical optimization task, the 
analyst must first define the objective along with any 
constraints to be imposed upon the objective variables; and 
finally a numerical algorithm must be prescribed to perform 
the optimization, preferably one which will do so effi- 
ciently for the particular problem. Since the term "optimum" 
most often refers to a solution obtained by mathematical 
programming, which is very inefficient for grid optimiza- 
tion, a near-optimum grid obtained by application of an 
adaptive procedure, henceforth will be termed an optimum 
grid. However, before such a grid can be determined to 
satisfy the stated objective of obtaining maximum solution 
accuracy for a prescribed cost, the terms "accuracy" and 
"cost" must be defined; but, more importantly, the optimiza- 
tion goals must be specified. This is critical because grid 
optimization can be implemented in various forms depending 
upon the optimization goals, which will, in general, be 
determined by the original purpose for performing the finite 
element analysis [Ref. 3]. For example, if the purpose of 
the analysis is to evaluate a local quantity, such as the 
maximum value of the dependent variable or one of its 



derivatives, then the nodal distribution should be densesr 
in the region where that taaximun is determined. If, on the 
other hand, the interest is on an integral quantity of the 
dependent variable ever a region of the domain, then the 
nodes should be assigned to achieve a uniform distribution 
of the error over that region. For the purpose of this 
investigation, attention will be focused on the three finite 
element resultants with the most significance in solid 

mechanics and other fields in which finite element analysis 

is employed; the maximum value of the dependent variable, or 
solution; the maximum value of the gradient of the solution; 
and the integral of the solution over the domain. 

In order to define the solution accuracy, it will be 
necessary to compare the error in the solution resultant 

obtained using an optimal grid to the error obtained using 
some baseline grid with the same number cf dsgrees-cf- 

freedem. For convenience, the reference grid chosen will be 
a uniform grid, or one with all elements of the same order 
and approximately the same size, with the understanding that 
no knowledgeable analyst would attempt to use such a grid in 
the solution of a problem with large gradients. 

The definition of analysis cost will be greatly simpli- 
fied by making the assumption that it is directly propor- 
tional to the number of d egrees-of- freedom used to obtain 
the solution. In reality the cost depends on many factors, 
some of which are very difficult to quantify. 

Understandably, the number of degrees-of-freedom is not the 
sole measure of computational costs, but it appears to be an 
adequate measure of preprocessing and postprocessing costs 
which often account for the major portion of the total 
analysis. 

This investigation will concentrate on the use of finite 
■element grid optimization methods for solving problems of 
structural mechanics. While this area has dominated the 
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literature on the subject, the techniques presented herein 



extend equally as well to any 
principles apply. 

There are two key questions 
to the adaptive application 
optimization : 

(1) Hhat criterion, based on 
finite element analysis, 
regions where the origina 



field for which variational 

which must be answered prior 
of finite element grid 

the results of an initial 
should be used to indicate 
grid is inadequate ? 



(2) What method of grid refinement should be employed ? 



Considerable attention will be devoted to these questions in 
the next two chapters. 

C. OBJECTIVES 

The objectives of this investigation are: 

(1) To examine some recently developed grid optimization 
techniques which have reached the state of practical 
application. 

(2) To implement some of these techniques in the solu- 
tion of some one- and two-dimensional linear self- 
adjoint boundary-value problems. 

(3) To draw a comparison among these applied techniques 
in terms of solution accuracy, analysis cost, and 
ease of implementation. 

(4) To recommend some guidelines to assist the analys* 
in obtaining optimal finite element solutions 
employing currently available or easily amendable 
software. 
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II. CRI TER IA FOR GRID RSFINE51M2 

The primary theoretical concern in the application of 
adaptive grid optimization is the selection of the refine- 
ment criterion. In ether words, one must decide upon which 
solution parameter, obtained from an initial idealization, 
may most appropriately be used to give some indication as to 
where the initial grid is inadequate and thus needs refine- 
ment. There are several competing proposals concerning the 
most appropriate choice of a refinement criterion. In 
reality, the decision must be based upon such factors as the 
type of problem being solved, which criterion is most prac- 
tically implemented, and whether accuracy is desired 
locally, globally, pcintwise, or with respect to an integral 
quantity. The following are some of the more practical 
refinement criteria used in grid optimization at present. 

A. SOLUTION PARAMETER 7AHIATION 

The most direct, computationally inexpensive, and 
earliest proposed indication of where an element grid 
requires refinement is a measure of the variation of some 
solution parameter over the domain. It is only logical that 
a piecewise polynomial approximation would experience the 
most difficulty in modeling the desired response in those 
regions where the solution or its resultants were either not 
smooth or were characterized by large gradients. Therefore, 
the basis cf this criterion is to refine the mesh in those 
areas where a solution parameter varies rapidly, with the 
implication that the optimum mesh is one for which the solu- 
tion parameter variation over each element is uniform 
throughout the domain. The first consideration in the 
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to find a schene for 



application of such a criterion is 
distributing the nodes to achieve such a condition. For 
one-dimensional problems the task is trivial, but one way to 
ensure equal variation over each elemeni; in higher dimen- 
sions is to position the nodes along equidistant contours of 
the chosen parameter. This is precisely the procedure 

prescribed for a practical optimization technique known as 
contouring. The other consideration is the determination of 
which solution parameter is to be used. In fact, several 
solution parameters have been found to work quite well 
[Ref. 4], but the most commonly used and the one that is 
consistent with the potential energy minimization formula- 
tion is the strain energy density [Ref. 2]. Because its 
employment requires only minor software changes and it has 
been found to produce excellent results, this refinement 
criterion was used extensively throughout the course of this 
research. 

B. GRID ITERATION 

Another, rather basic but less direct, method of 
locating regions of the mesh to be refined is known as grid 
iteration, which can be implemented in one of two ways. An 
initial coarse grid analysis may be repeated using either a 
finer or a higher order mesh, and a comparison of the resul- 
tants of interest between the two solutions will identify 
those areas of the domain where refinement is most effec- 
tive. Another approach is based on the assumption that the 
greatest benefit is to be gained by refinement in those 
regions where a small perturbation, like the introduction of 
a single additional degree- of-freedom , produces the greatest 
change in the solution or one of its resultant parameters. 
Since additional degrees-of-f reedom would be expected to 
produce the greatest change in those regions where the 
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desired response varied most rapidly, refinements based on 
this method provide results very similar ro those obtained 
using the solution parameter variation criterion already 
discussed. The grid iteration method may at first appear to 
be more computationally expensive, but it was developed to 
be most efficiently implemented employing a special family 
of elements. These elements, called "hierarchical", possess 
some very desirable properties for this application and will 
be discussed in the next chapter. 

C. EL3HENT RESIDOALS 

The major drawback with refinement criteria based on 
solution parameter variations is that their applicability 
appears limited to elastostatic problems. For this reason, 
several investigators have recently developed grid refine- 
ment criteria based on element residuals, which appear prom- 
ising for application to all types of finite element 
problems, including nonlinear analysis. The reason for this 
is that the residual has essentially the same meaning 
regardless of the problem type rsef. 5]. For example, 
consider the governing differential equation, 

D C u ] - f = 0 

defined on some domain, where D [ ] is a linear or nonlinear 

differential operator, and the dependent variable u and the 
non- homogeneous term f are both functions of the independent 
variables. Let the finite element approximation to the 
solution of the differential equation be u - u. Applying the 
differential operator to the approximation gives rise to the 
residual, which is defined as 

R = D [ U ] - i 

and is not identically zero unless the finite element solu- 
tion is exact. 
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The key to using the residual as a criterion for deter- 
mining regions of the domain where refinement is necessary 
is the local residual on the element level, which indicates 
the contribution of the element to the rotal error of the 
finite element approximation. Since the residual is a point- 
wise quantity, the useful measure of the element error 
contribution is a norm of the element residual, or the inte- 
gral over the element of the product of the residual and 
some weight function. The integration is most readily 
performed using numerical quadrature. The grid optimization 
strategy then becomes one of refining the mesh so as to 
egui -distribute the element residual norms, by forcing them 
to become smaller and more uniform over the domain. 

There are some drawbacks to the element residual refine- 
ment criterion which have not yet been fully resolved, such 
as appropriate selections of the residual norm and the 
weight function, and in the computation of rhe residual 
itself. While the evaluation of the residual presents no 
particular difficulties in the interior of the element, iz 
is rarely determinable at the edges. The reason for this is 
than in formulating the finite element approximation the 
element shape functions are defined so as zo provide only 
that degree of continuity required to adequately model the 
physical problem; the most frequent consequence being that 

'\j 

D [u] is undefined along the interelement boundaries. 
Unfortunately, this singularity cannot be ignored and a more 
complicated analysis must be applied in order to bound the 
residual contributions at these boundaries [Ref. 6]. 

D. A-POSTEBIOai ERROR ESTIMATES 

A sophisticated extension of the element residual 
criterion is one based on computable error estimates from an 
initial finite element analysis. This utilizes the energy 
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norm of the residual, in which case the weight function is 
the residual itself. Research in reliable error estisiates 
was pioneered by Babuska [Ref. 7, 3] for linear quadrilat- 
eral elements and more recently by Zienkiewicz [Ref. 9] for 
higher order elements. The major difference from the 
residual criterion is that instead of egui-distr ibuting the 
element residual norms over the domain, they are normalized 
to compute error indicators for the elements, which are in 
turn used zo compute reliable poiatwise error estimates for 
the solution as well as the energy error over the domain. 
These quantities are of primary importance because rhey 
provide not oaly an indication of where additional refine- 
ment is most effective, but also a measure of the quality of 
the mesh no determine whether additional refinement is 
necessary [Ref. 9]. The optimization strategy is to obtain 
a nearly uniform distribution of the error indicators 
throughout the domain, which corresponds to minimizing the 
error in the energy norm. The refinement procedure may prog- 
ress unril all the local errors are within some prespecified 
tolerance. While the practical utility of such a refinement 
criterion is obvious, the marhematical developmenr and the 
algorithms involved are rather complicated. However, the 
process is not compute rionally expensive, and there now 
exists a prototype self-adaptive finite element code, FEARS, 
which implements this refinement criterion [Ref. 10]. 



18 



III. METBODS OF GRID REFINEMENT 

Once it has been determined where the initial element 
grid is inadequate and needs refinement, the next considera- 
tion is how the idealization in these areas should be 
improved. The choice of the refinement method to be employed 
may well be a more important decision than the selection of 
one of the refinement criteria previously discussed, since 
at least one investigator has observed that for a particular 
method of grid refinement, the various refinement criteria 
produce essentially the same solution results [Ref. 11]. 

Grid refinement is the process of introducing additional 
degr ees-o f- freedom into selected regions of the finite 
element grid, and may be performed by one of three methods; 

(1) The polynomial degree of the elements remains fixed, 
usually at a low value, while the size of the 
elements is reduced. This has become Icnown as the 
h-version since element size is commonly denoted by 
the letter h. 

(2) The size of the elements, usually few in number, 
remains fixed while the polynomial degree of the 
elements is increased. This has become hnown as the 
p-version since polynomial order is commonly denoted 
by the let*er p. 

(3) The size of the elements may be reduced concurrently 
with an increase in their polynomial order. This is 
known as the combined h- and p-version of the finite 
element method. 
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A. C0HTEBG2NCE OF GBID REFINEMENT 



It is well known that the finite element method 
converges with an increasing number of degr ees-of- freedcra; 
in fact, this is the justification for its development. 
Therefore, the appropriate measure of the effectiveness of a 
particular grid refinement method should be the associated 
rate of convergence, which generally will be affected by the 
smoothness of the approximated function over the subdomain 
of interest. It has been demonstrated that when the refine- 
ment is performed uniformly over the domain, the associated 
rate of convergence is asymptotic, provided the number of 
degr ees-of-freedom is sufficiently large [Ref. 1]. The 
asymptotic rate of convergence is often measured as the 
slope of the error versus cost curve in the linear, or 
asymptotic, range when plotted on a log-log scale. In 
performing such an error analysis for the displacement 
formulaticn of the finite element method, the error is 
usually the relative strain energy error, approximated by 
the energy norm, and the cost is assumed to be some simple 
function of the average element size or the number of 
degrees-cf- freedom [Ref. 12; p. 726]. Only in the past 
several years has there been any significant research 
comparing the relative merits of the different methods of 
grid refinement. Since the solutions of elliptic boundary 
value problems are usually very smooth over convex domains 
except in the vicinity of corners, most of this research has 
focused on solutions exhibiting singularities, which 
severely hinder the rats of convergence, as in problems of 
fracture mechanics and in domains with re-entrant corners 
[Ref. 1 , 13, 14, 15]. 

In order for a finite element analysis to be both effi- 
cient and reliable, the asymptotic convergence range should 
be entered for as few degrees-of -freedom as reasonably 
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possible. In general the p-version satisfies this require- 
ment better than the h-version. While it has been estab- 
lished that p-ccnyergence will nsoessarily occur whenever 
h-convergence occurs, the converse is not true. For example, 
if the h-version using a uniform grid of linear elements is 
applied to a nearly degenerate problem, the number of 
degrees- of- freed cm required for entry into the asymptotic 
range may be beyond the computer's round-off limitarions, in 
which case caivergence will not occur unless rhe polynomial 
order is increased [Ref. 1]. Numerical experiments on such 
problems clearly indicate that the p-version requires 
considerably fewer degrees- of-freedom than the h-versicn to 
achieve the same degree of accuracy. Recent analyses 
[Ref. 1, 13] of the asymptotic rats of convergence in energy 
for non-smooth solutions, using uniform refinemenz with 
sufficiently high numbers of degress -of -freedom, have demon- 
strated that the p-version cannot have a slower rate of 
convergence than the h-version. Furthermore, if the singu- 
larity is confined to element boundaries, as is usually the 
case, the error for p-raethod is inversely proportional to 
the number of degr ees-of- freedom, whereas the error is 
inversely proportional to the square root of the number of 
degrees-of-freedom in the h-version. In other words, for 
this special class of problem, the rate of convergence for 
the p-version is twice that for the h-version, which is due 
primarily to the ability of higher order polynomials to 
"absorb" singularities occurring at the element boundaries. 
This implies, at least for this type of problem, that in 
order to minimize the error for a specified number of 
degrees-of-freedom, the best strategy is not to subdivide 
the domain uniformly, but to use instead a single element of 
increasing polynomial order [Ref, 15]. 
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since it is unlikely that ona would antsiapt to solve 
such a problem using uniformly finer grids, a more useful 
comparison between the convergence rates of the two versions 
would be based on adaptive refinement employing one of the 
solution -based criteria discussed in the previous chapter. 
It so happens that the h-version, when used with optimally 
refined meshes, can have a higher convergence rate than the 
uniformly distributed p-version, provided that the element 
order is sufficiently high. However, the p-version can also 
be employed with an optimal refinement criterion. While 
there are yet no proven theorems concerning the convergence 
rates for non-uniform refinement, obtaining the desired 

solution accuracy with optimal p-dis tributions appears to be 
much less sensitive to the particular choice of the elements 
to be refined than with optimal h-refinement [Sef. 13], 

It would seem plausible that an even better optimization 
strategy would involve a proper combination of both the h- 
and p-versicns. It has been demonstrated for problems with 
corner singularities, that a proper sequence of 

h-ref inements combined concurrently with the proper sequence 
of p-distributions results in extremely high convergence 
rates, conjectured tc be exponential [Sef. 15], However, 
this proper combination is difficult to determine, and adap- 
tive refinement based on the combined h- and p-versions 
poses some difficult data management problems. To avoid this 
problem a more promising approach, proposed by Babuska and 
Szabo [Ref. 1], employs a graded mesh in which the element 
sizes are first reduced according to a geometric progression 
towards the singularity, followed by determining the optimal 
p-distribution for those elements using an adaptive 
criterion. However, obtaining the optimal combination when 
employing this scheme can be a delicate matter and, astcund- 
ingly, the highest accuracy is achieved when the polynomial 
order of the elements actually increases with distance from 
the singularity. 
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There are sciae additional advantages of the p-version 
worth mentioning. Because the p-version employs fewer 
elements, there are lesser preprocessing and postprocessing 
costs than for the h-version. Furthermore, when bandwidth 
minimization and sparse matrix solution techniques are used, 
the solution time for the p-version is approximately the 
same as for the h-version for a specified number of 
degrees- of- freedom, and the p-version appears less suscep- 
tible to round-off errors. Finally, the p-version is simpler 
to implement adaptively than the h-version when hierarchical 
elements are employed [Ref, 13], 



B. HIERARCHICAL FINITE ELEMENTS 

The hierarchical concept was first introduced as a 
simple method for implementing the p-version and as a 
convenient device for imposing boundary continui-y between 
elements of different polynomial order [Ref. 9], Since then 
a useful family of elements based on the hierarchical 
concept has been developed and incorporated into COMET- X, an 
experimental finite element code, developed by Szabo, which 
self -adaptively employs both the h- and p-versions of the 
finite element method [Ref. 14], 

For a brief description of the hierarchical concept 
consider the conventional finite element formulation which 
produces the following system of equations: 



^(n)- f(n) 



K, ,u^“^= f\ 



(Eqn 



. 1 ) 



where n is the number of degrees-of- freedom, is the 

n X n global stiffness matrix, u^^^ is the finite element 






approximation of the exact solution, 
n-component global load vector. When 



and 



f(n) 
hiaher 



the 
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dsgr e€s-of-frs6dcm are added to the original system using 
conventional refinement methods nhe system of equations 
becomes: 




f (n-tai) 



(Egr. 3.2) 



where the contributions to K 



(n-fm) 



and f^’^'*™^from the refined 






elements result in a completely different set of coeffi- 
cients. If, on the other hand, this refinement had been made 
hierarchically, the equations would become: 
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f(n) 


f^(n) f^(n,m) 
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where and f^^^jare the stiffness matrix and force vector 

from the original system of equations for n degrees-of- 
freedom appearing in Equation 3.1. However, u^^^is not the 
nodal values of the finite element solution for the m 
additional degrees-of-freedom, but instead represents the 
difference between those values and the pointwise values 
obtained from the lower order polynomial interpolation for 
the unrefined mesh of n degrees-of-freedom. 

The primary advantage of hierarchical elements is imme- 
diately observable from Equation 3.3. Because the shape 
functions of an element of order p consxitute a subset of 

the shape functions of an element of order p+1, the local 
stiffness matrix and force vector for each hierarchical 
element is embedded in the stiffness matrices and force 
vectors of all hierarchical elements of higher order. 
Therefore, the global stiffness matrix ^^^^and force vector 
f^'^^ of the original system are preserved, thus saving 
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considerable time and effort expended on computing the coef- 
ficients for successive refinements [Ref. 14], Another 
advantage is that the hierarchical form of the global stiff- 
ness matrix is more diagonally dominant than the one 
resulting from a conventional refinement, resulting in 
improved conditioning and faster convergence when iterative 
solvers are employed [Ref. 9]. Another benefit of hierar- 
chical elements, which arises from the "add-on" nature of 
the nodal variables of the higher order degrees-of-freedcm, 
is that the problem of maintaining boundary continuity 
between elements of different polynomial order becomes 
trivial. Instead of introducing global constraint equations 
for the higher order degrees-cf-freedom, the nodal variables 
are simply set equal to zero and condensed out, as if they 
were zero-valued Dirichlet boundary conditions [Ref. 2]. 

There are two major drawbacks with hierarchical elements 
that have hampered their widespread acceptability. The 
first, which has already been mentioned, is that the nodal 
variables for the higher order degrees-cf-freedom represent 
difference values rather than the more easily identifiable 
values of the dependent variable itself. Secondly, when 
implementing the h-version of the finite element method, 
special integration rules must be introduced when the subdi- 
vided element is in hierarchical form [Ref. 9]. Of course, 
the latter problem can be evaded by using the p-version, fcr 
which the hierarchical concept was developed. In spite of 
the disadvantages of hierarchical elements, their consider- 
able computational efficiency and utility for grid optimiza- 
tion will certainly result in their widespread utilization 
in future adaptive finite element codes. 
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IV. GRID OPTIMIZATION TBCHHI^OES 



Once the analyst has identified where the initial grid 
needs enrichment and decided which refinement method to 
employ, he must then determine a sysremaric procedure, or 
algorithm, to perform the refinemeni. according to the 
criterion selected. The ultimate goal of such a procedure is 
to design an element grid which meets the optimization 
objective of obtaining maximum solution accuracy for a spec- 
ified analysis cost. While the analyst may or may not have 
an indication of the accuracy of the solution, he should 
have a preconceived notion of cost, or how much effort he is 
willing to expend to arrive at a better solution. Therefore, 
with some Icnowledge of the grid optimization techniques 
available and an understanding of the advantages and disad- 
vantages of each, the analyst can realize the grid optimiza- 
tion goal. 

There are essentially two adaptive grid optimization 
strategies: 

(1) Grid refinement, in which the initial analysis is 
performed on a relatively coarse grid, and new 
degrees-of -freedom are added to the same grid by the 
iterative application of the solution-based 
criterion. 

(2) Grid modification, in which the initial analysis is 

performed using a prsspecified number of 

degrees-of-freedom, and the solution-based criterion 
is employed to shift degrees-of-freedom from certain 
regions to others. This may involve complete grid 
redefinition in an effort to obtain a near-optimum 
grid in a single cycle. 
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Much of the interest lately has been in the develop^ient 
cf complicated self-adaptive software packages which mini- 
mize the impact of the user’s skill on rhe final solution. 
Ideally, the analyst would merely define the problem and the 
program would automatically generate and analyze the optimum 
grid employing one or more of these techniques, possibly at 
the user’s option. 



a. MATHEMATICAL PROGBAMMIMG 

No discussion of grid optimization techniques would be 
complete without a brief description of mathematical 
programming, not only because it is how grid optimization 
was earliest attempted, but more importantly, it is 
precisely what the engineer envisions when ha hears the term 
"optimization". It is nor a grid optimization rschnigue, per 
se, but rather a numerical process of achieving any optimi- 
zation objective, once it is explicitly defined in mathemat- 
ical terms. In solid mechanics the finite element method is 
a numerical method for minimizing the potential energy func- 
tional, which in discretized form may be written; 



IT = k u''' K u 



f 

a, r\j 



(Egn. 4.1) 



where : 



u is the global displacement vector 

K is the global stiffness matrix, and 

% 

f is the global load vector 



In the classical finite element 
energy is minimized with respect 
which implies satisfaction of 



formulation, the potential 
to the nodal displacements, 
the following stationary 



conditions: 
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0 



.^n) 



(Egn. 4.2) 






3ui 



(i = 



where n is the number of degr eas-oi- freedom. This leads to 
the very familiar system of linear aquations: 



^ a - f =0 (Eqn. 4.3) 

However, since ^ and f ar e functions of the nodal coordi- 
nates, then the potential energy could be minimized with 
respect to the nodal coordinates as well. This would require 
satisfaction of the following additional stationary 
conditions: 

3TT 

3xj = 0 (j=1,2,...,a) (Sqn.4.4) 



where m is the number of nodal coordinates, . This differ- 
entiation leads to the less familiar system of non-linear 
equations: 



T 

'1 2^’ u - u = 0 

Oi 3 xj c»Xj % 



(j = lf2,. 



i m) 



(Egn. 4.5) 



This, then, is the mathematical statement of the grid opti- 
mization problem for the elastostatic case. The nodal 
displacement, variables may be eliminated by minimizing the 
potential energy with respect to the nodal coordinates only, 
subject to the implicit constraint that Equation 4.3 is 
always satisfied [Ref. 4], Unfortunately this does not help 
much because the objective function is still nonlinear, 
rendering most numerical optimization algorithms inefficient 
and unreliable. The difficulty is even further compounded by 
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the requirement that the nodal variables be subject to side 
constraints in order to maintain the defined boundary of the 
domain and to ensure that the elements neither distort 
excessively nor overlap one another. For all except the 
simplest of problems, these constraints may be even more 
severely nonlinear than the objective function, resulting in 
the analysis becoming prohibitively expensive [Ref. 2]. For 
this reason, mathematical programming in finite element grid 
optimization has been abandoned in favor of some equally 
reliable, yet far mere computationally efficient grid opti- 
mization techniques. However, these early efforts with 
mathematical programming were nor tonally in vain because 
they gave rise to the contouring techniques. 

B. COHTOOBINS 

Since mathematical programming is infeasible for grid 
optimization, further in vestigari ons were conducted to 
suggest some guidelines to enable rhe analysr to construct a 
grid with similar topological features of the numerically 
optimized grid without the computational effort involved. 
Turcke [Hef. 4], in employing mathematical programming in 
the solution of some simple two-dimensional elastcstatic 
problems, observed xhat there was a very definire element 
pattern common among problems involving high strain gradi- 
ents and that the nodes of the numerically optimized grid 
generally tended to be aligned along contours of some 
response function being modeled. Consequently, in performing 
analyses on grids whose construction was based on contours 
derived from an initial analysis, it was determined that the 
following provided grid characteristics in regions of high 
strain gradients similar to the numerically optimized grid, 
but at a fraction of the computational expense: 
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• contours of displacement 

• contours of maximum principal stress 

• contours of maximum shear stress 

• contours of strain energy density (isoenergetics) 

• principal stress trajectories (isostatics) 

Since the strain energy density is the response which is 
consistent with the principle of minimum potential energy, 
isoenergetics are the most commonly used contours along 
which element edges are aligned [Baf. 4 ]. However, there 
still remains the question of how to position the nodes 
along the contours. For this reason, isostatics have become 
increasingly popular because the principal stress trajecto- 
ries form a "flow net" of orthogonal curves which can guide 
the analyst in the layout of the elements [Ref. 16]- 

Since contouring involves the redefinition of the grid, 
as opposed to a grid refinement, its primary advantage is 
that the enriched mesh is not constrained to the element 
configuration of the previous mesh. Therefore, there is no 
limit to the amount of enrichment per cycle which can be 
performed and it is conceivable that an optimum mesh could 
be generated in a single cycle [Ref. 2]. However, while the 
computational cost of repeated analyses is reduced, the 

preprocessing costs involved in constructing the contours 
and redefining the mesh can become quite high, especially if 
the contours are complex. Algorithms for performing these 
tasks in two-dimensional domains have been proposed 
[Ref. 4, 11], but they are not extendable to three- 
dimensional problems. The major obstacle for two- and 
three-dimensional domains is that it is often difficult to 
constrain the element edges to the contours without the 
elements becoming elongated or distorted to the degree that 
numerical inaccuracies result. Another difficulty, not 

addressed in the literature, is how the contour increments 
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should be selected when the response function is non- 
monotonic over the domain. 



C. SELECTIVE REFINEMENT 

The most commonly employed grid optimization rechnigue 
is that of selective refinement. As its name implies, 
selected elements from a given mesh are enriched while the 
original element grid remains essentially intact. The 
elements selected for refinement are determined by the iter- 
ative application of the solution-based criterion to indi- 
cate which elements contribute most to the solution error. 
The refinement can be performed by either the h-versicn or 
the p-version, or even the combined version if so desired, 
but the choice is most often predetermined by the capabili- 
ties of the available preprocessor. Since the addition of 
new degrees-of -freedom over several iterations can guickly 
enlarge the problem, it is advisable to perform the initial 
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Figure 4.1 Seme h-7ersion Subdivision Schemes. 

the case in the h-version scheme of Figure 4.1 and ir also 
arises in the p- version when two elemenrs of different poly- 
nomial order share a common edge. when rhis situation 
occurs, the additional de grees-of- freedom do not actually 
represent degrees-of-freed om an all because they must be 
numerically constrained to the polynomial interpolant of the 
lower order. Such censrraints are usually imposed in one of 
three ways: global constraint equarions may be writren; the 
constraints may be incorporated in the elemental basis; or 
hierarchical forms may be used wirh the excess degrees-of- 
freedem simply set to zero and condensed out [Ref. 2]. 

There are some other selective refinement techniques 
which do not require any major software modifications. In 
the h-version, the continuity problem may be circumvented by 
employing any of the coarse-to-fine mesh transition schemes 
for which all of the element edges remain of the same poly- 
nomial order [Ref. 17: p. 210]. However, it is impossible 
to employ these schemes without permitting some element 
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distortion, and the refinement must nearly always be 
performed interactively rather than automatically. For the 
p-version, interelement continuity can be easily ensured by 
employing vacia ble-noded isoparametric elements, which 
permit a single element to possess edges of differing poly- 
nomial orders [Ref. 17: p, 125]. 

The analyst must also exercise care when adding new 
nodes to the boundary of the domain to ensure that the 
appropriate boundary conditions are determined and applied. 
Furthermore, if the boundary is curved, the coordinates of 
the new node should be computed such that it is placed on 
the actual boundary and not necessarily on the edge of the 
element being refined [Ref. 2]. 

The important advantage of the selective refinement 
technique is that once an appropriate refinement criterion 
has been determined, selecting candidate elements for 
refinement in each cycle becomes straightforward. The 
refinement can then be continued indefinitely to achieve 
very high accuracy, but because the solution phase is 
repeated for each cycle, it is desirable to hold the number 
of cycles to a minimum. Because the nodes from the previous 
mesh remain fixed for each cycle, selective refinement is 
ideally suited for iterative solution methods. The solution 
values obtained from the previous cycle, combined with 
interpolated values for the new degrees-of-freedom, provide 
an excellent initial guess for the next cycle [Ref. 2]. 

The major disadvantage is that the limited amount of 
refinement which can be performed in each cycle may necessi- 
tate several cycles to obtain an optimum grid. In addition, 
if new degrees-of-freedom require interelement continuity 
constraints, data management can become cumbersome unless 
the constraint is performed hierarchically. 
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D. SOBDOHAIN ISOLATION 



One of the obvious disadvantages of the selective 
refinement technique is that the solution must be completely 
repeated for each cycle when, in fact, rhe number of 
degrees-cf- freedom added in each cycle may be few in compar- 
ison to the total for the problem. In addition, the number 
of elements requiring refinement in each cycle may only 
account for a small portion of the domain. Although the 
refinement criterion has indicated where the grid is inade- 
quate and the approximation is likely to be poor, the solu- 
tion is repeated in each cycle for those nodes where the 
error is presumably small. Besides the apparent computa- 
tional inefficiencies, this shortcoming severely restricts 
the amount of refinement which can be performed in the 
subregions of interest since it is desirable to confine the 
size of the problem within reasonable limits. An alxernative 
approach is to reformulate the problem for those subregions 
where refinement is necessary and to accept the results of 
the initial analysis as an adequate solution for the 
remainder of the domain. The elements requiring further 
refinement, which constitute isolated subdomains of the 
original problem, can generally be subjected to signifi- 
cantly greater refinement than would otherwise be practical. 
The solution obtained from the initial analysis is then used 
in imposing boundary values on those degrees -o f-freedcm 
located on the boundaries of an individual subdomain. These 
can, in turn, be used to generate the boundary conditions 
for any additional boundary degrees-of-freedom introduced by 
the refinement using an appropriate interpolation scheme. 

This grid optimization technique, which the author terms 
"subdomain isolation", has some further advantages over 
selective refinement. The subdomain may be selected arbi- 
trarily small such that excellent results may be obtained 
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with, a single cycle using uniform refinement. Therefore, “he 
difficulties involving coa rse-to-f ine transition schemes, 
element elongation and interelement continuity can be 
avoided. Furthermore, one can choose as many subregions for 
refinement as desired without creating an excessively large 
problem. 

The obstacle which may prevent this technique from being 
readily accepted is the notion that, by imposing erroneous 
boundary conditions on the subdomains, the convergence of 
the finite element method to the exact solution in these 
regions has somehow been tampered with. This aversion may be 
somewhat abated by considering a simple extension of 
Saint-Venan t 's Principle [Ref. 18: p. 33], Although the 

conditions are not rigorously satisfied at the boundary, 
which may result in significant changes in the response 
locally, the effect at some sufficient distance away will be 
negligible. The numerical evidence supports this premise. 
While errors in the boundary values may somewhat restrict 
the accuracy cf the dependent variable, great improvements 
can be realized in the accuracy of its gradients, which is 
more often the goal of the optimization. Since the initial 
analysis provides the boundary values for the subdomains, it 
is desirable that its solution be as accurate as reasonably 
possible. Fortunately, since subsequent refinements are not 
performed on the original grid, the initial analysis may 
involve significantly more degrees-of-freedom than in the 
case of selective refinement. 

E. SESH GRADING 

The final grid optimization technique to be discussed 
employs a mesh for which the element sizes are successively 
reduced, according to some geometric sequence, towards a 
selected region of the domain. One might argue that mesh 
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grading is not really an o ptimizacio n technique since it is 
most often applied on an ” a-priori" basis rather than adap- 
tively, and that it does not lend itself well to "he itera- 
tive application of a solution-based criterion. However, the 
technique is simple to use and its implementation requires 
few software modifications. Furthermore, a solution-based 
refinement criterion can be used to give a measure of the 
quality of the mesh to indicate waether a more pronounced 
grading may prove beneficial. Depending on the solution 
parameter of interest, mesh grading can provide excellent 
accuracy at a low analysis cost. This refinement method must 
therefore be considered among the grid optimization 
techniques. 

For the less elaborate finite element preprocessors, 
mesh grading is often the only refinement means available 
without resorting to a uniformly finer mesh involving many 
more degrees-of- freedcm. The most common method of implemen- 
tation in two-dimensions is to first define the problem 
domain in terms of a curvilinear quadrilateral by selecting 
four keynodes along the problem boundary. Then the boundary 
nodes are spaced according to some geometric sequence based 
on the user-provided bias parameters which determine the 
density of the nodes towards selected points on the four 
quadrilateral edges. Finally, carves are generated to 
connect the boundary nodes on opposite edges of the quadri- 
lateral, thus producing a graded mesh. This process, which 
can also be extended to three-dimensions, is the mesh gener- 
ation scheme employed in the finite element code GIFTS 
[Ref. 19]. 

The major disadvantage of mesh grading is that in order 
to achieve sufficiently small elements in the region of 
interest, the elements must grow successively larger away 
from that region. This may be very undesirable, especially 
if refinement is called for in more chan one region of the 
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domain, in which case the mesh must be 
by subdomains, thereby complicating 
involved. 

Another disadvantage is that unless 
some special geometric symmetry. 



generated and graded 
the data management 

the domain possesses 
excessive element 




Figure 4.2 Graded aesh for a Perforated Square Plats. 

elongation will usually result if a highly pronounced 
grading is required, some configurations are particularly 
well suited for refinement for mash grading such as the 
classical perforated square plate problem in solid mechanics 
shown in Figure 4,2. 
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APPLICATION TO ONE-DIMENSIONAL PROBLEMS 

Now that the necessary tools for performing grid optimi- 
zation have been introduced, it is tine to employ them in an 
attempt to obtain optimal solutions to some practical prob- 
lems in engineering. An obvious starting point for such an 
investigation is the one-dimensional boundary value problem. 
While most of the fruitful research in grid optimization has 
concentrated on problems of higher dimensions, the one- 
dimensional problem is a very convenient device for studying 
finite element grid optimization. Foremost, one-dimensional 
finite element models possess a unique connectivity in that 
adjacent elements meet at their end nodal points. Therefore, 
refinement by the h- or ?-versions, or by relocating nodal 
points becomes a trivial task, which does not involve any of 
the difficulties so frequently encountered with higher 
dimensional problems, such as preserving interelement conti- 
nuity and maintaining optimal element shapes. Furthermore, 
one- dimensional studies can often provide valuable insight 
to the solution of more difficult higher dimensional 
problems. 

The primary concerns in the selection of the problems to 
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Because of the complexity and a certain degree of arbi- 
trariness involved in the computation of element residuals 
and a-posteriori error estimates, the solution parameter 
variation is the refinement criterion of choice. There are 
several solution parameters which are easily computed, 
requiring miaimal software changes to an exisring finite 
element code. 

Furthermore, for the one-dimensional investigation, ir 
was decided to simplify the analysis by exploiting the 
linear elements. While it is granted that improved solution 
accuracy may generally be obtained by employing higher order 
elements, it will be assumed that conclusions based on the 
use of linear elements can be applied as well to elements of 
higher pclynomial order. 

A. ELASTIC CABLE PROBLEM 

Consider an elastic cable under tension r, stretched 
between two points a distance 2L apart. If the cable is 
supported by a Winkler, or elastic, foundation of modulus k, 
and a concentrated load P is applied at the midspan, the 
resulting deflection v(x), (0 < x < L) , is as shown in 

Figure 5.1. The analytical solution and finite element 
approach for this problem are presented in Appendix A. 

The initial finite element analysis was performed using 
ten linear elements of uniform length. From this initial 
analysis the approximate distribution of the following solu- 
tion resultants was obtained over the domain: 

• the displacement, v (the solution) 

• the slope, v' 

• the strain energy, U 

• the strain energy density, SED (dU/dx) 
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Figure 5.1 Tension Cable Deflection on Elastic Foundation. 

Subsequent analyses were performed for finite element models 
using the same number of elements, but with the nodes redis- 
tributed to achieve approximately uniform variation of the 
above parameters over each element. Note that the strain 
energy refinement criterion produces elements of identical 
strain energy content. In addition, the problem was solved 
employing graded element models of various adjacent element 
length ratios. The resulting element models based on these 
refinement criteria are shown in Figure 5.2 (a-f). The 

graded model (b) fcr an element length ratio of 1.2 is 
presented for comparison because it produced good overall 
solution results. 

As previously mentioned, the solution resultants of 
primary interest are the maximum displacement, the maximum 
slope, and the integral of the displacement over the domain, 
because they represent important analogous solution results 
in nearly all fields in which finite element analysis is 
often performed. The accuracy obtained in these values for 
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(a) Uniform 



(b) Graded 1.2 



(c) V 




(d) V* 



(e) U 



(f) SED 



^mod 




(h) 



SED 



mod 



Figure 5.2 Tension Cable Refineaents 




each of the refined models is presented in Table I. As can 
be seen in Figure 5,2, the strain energy and strain energy 
density criteria produced extreme variations of element 
length while the criteria of displacement and slope result 



TABLE I 

Tension Cable Problem Solution Results 



Problem Parameters; L = 1 00 in 

k = 1 psi 

T = 1000 lb 
P = 1000 lb 





Variation 

Refinement 

Criterion 


Percentage Relat 
V (max) V ' (max) 


ive Error 
[vdx 




Uniform 


-0. 40 


0. 36 


0. 12 




Graded (1.2) 


-0. 19 


0. 07 


0.17 




V 


-0. 18 


0. 06 


0.39 




V' 


-0.23 


0. 05 


0.87 j 




U 


-1 . 03 


0. 05 


3.58 




SED 


-1.29 


0.05 


4. 17 




U (mod) 


-0. 53 


o.ou 
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SED (mod) 


-0.51 


0.03 


1 .48 

j 
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therefore refinement cannot reduce its error. fst great 
improvement in the accuracy of the maximum slope and modest 
improvement in the accuracy of the maximum displacement can 
be achieved with moderate refinements based on the displace- 
ment and slope distributions. 

One might assume, and correctly so, that the ability of 
the energy refinements to produce the best accuracy for the 
maximum slope is due to the extremely small elements which 
result in the area where that quantity occurs. Furthermore, 
it would be correct to propose that the reason for these 
refinements producing poorer estimates than the uniform 
model for the other two quantities of interest is that the 
excessively large elements in the regions of low gradients 
severely overstiffen the model there. It would then seem 
plausible tc improve the accuracy for the maximum displace- 
ment and the integral quantity by redistributing the nodes 
in these regions to prevent such excessively large elements. 
This was done by arbitrarily employing a grading scheme to 
the three largest elements to produce the modified refine- 
ments based on strain energy and strain energy density shown 
in Figure 5.2 (g) and (h) . As can be seen in Table I, such a 
modification did indeed significantly reduce the errors in 
the maximum displacement and the inregral, but it even 
further improved the accuracy for the maximum slope. 

One might conclude from Table I that the best overall 
model was obtained using the graded mesh, and that since it 
is easier to obtain, it should be deemed the optimal grid. 
But this particular grading was cnosen for precisely chat 
reason and was presented only as a means of comparison. In 
practice, the selection of a grading ratio is somewhat arbi- 
trary and making an adequate choice may be difficult. 

There is justifiable confusion as to which refinement 
produced the "best” solution accuracy for this problem and 
it raises perhaps the most important issue in the subject of 
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grid refinement. Before any optimization process can be 
pursued, the optimization goals must be explicitly defined- 
Clearly, as is the case in this problem, the designation of 
the optimum grid would depend heavily upon which of the 
three solution resultants is most critical to the analysis. 

B. TiPEBED BAB PBOBLEB 

The linearly tapered bar under axial loading has 
received considerable attention and was one of the early 
problems for which analytical grid optimization was 
employed. Consider a tapered elasric bar of length L and 
modulus E, fixed at one end, with an axial load P applied at 
the other, . for which the axial displacement u (x) , 
(0 < X < L) , is desired. The cross-sectional area varies 
linearly from at the fixed end to at the tip, as shown 




Figure 5.3 Linearly Tapered Bar Under Axial Loading. 

in Figure 5.3. The analytical solution and finiae elemena 
approach are presented in Appendix 3. 
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One of the significant features of the tapered bar 
problem is that the maximum stress can be very difficult to 
model accurately, and it is for precisely these problems 
exhibiting large strain gradients that grid optimization 
becomes most beneficial. Interestingly, the stresses 
obtained at the element midpoints are exact for this 
problem, and the difficulty arises from the inability of the 
constant slope shape functions to model the maximum stress 
occurring at the boundary. In examining this problem, 
Prager [Sef. 20] demonstrated analytically that when each 
element has the same strain energy content, the relative 
error in displacement is identical for all the nodes. 
However, this phenomenon appears peculiar to this problem 
and the author does not subscribe to such a measure of an 
optimum grid. Judging the effectiveness of a particular 
refinement based upon the deviation or the mean value of the 
pointwise errors generally tends to be unfavorable to opti- 
mization procedures since they almost always introduce many 
more nodes in those regions where the response is most 
difficult tc models Hence, an improved solution may have a 
larger mean value of the pointwise errors [Ref. 3 ]- 

The criteria employed in the refinement of the tapered 
bar model are identical to those used in the cable problem 
and their effects are shown in Figure 5,4 (a-e) . Two excep- 
tions are that now the displacement and strain energy 
criteria produce identical refinements, and the graded model 
chosen as the best overall is now based on a grading ratio 
of 1.4, producing a more drastic refinement than that of 1.2 
for the cable. This further demonstrates the difficulty 
involved in obtaining adequate element grids on an 
"a-priori" basis. 

The solution results are presented in Table II and the 
most readily apparent observation for this problem is the 
large errors in the maximum slope, which would severely 
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(a) Uniform 








(f) 



SED ^ 
mod 



Figure 5.4 



j.apered Bar Sefinemen'cs. 



46 



TABLE II 

Tapered Bar Problem Solation Results 
Problem Parameters: L = 100 in 





A = 13.5 in2 
A = 0.5 in2 
E = 10x10* Dsi 
P = 10x103 lb 






Variation 

Refinement 

Criterion 


Percentage Relative 
u(max) u* (max) 


Error 




Uniform 


-3.80 


-37.5 


0.68 




Graded (1.4) 


-0.78 


-4 . 1 


0.14 




u ; U 


-0. 85 


-10.6 


0. 15 




u* 


-1.81 


-7.7 


0.33 




SED 


-6. 54 


-3.6 


1.18 




SED (mod) 


-1.99 


-3.6 


0. 36 


■ 


underestimate the maximum stress. 


These resul 


ts are 


based on 


quadratic extrapolation 


of the exact slopes 


at the 


element 


midpoints, since the 1 


inear shap 


9 functions 


would 


produce 


even poorer estimations 


of the max 


imum slope. 


As before, the 


more extreme refinement 


based on 


the strain 


energy 


dans ity 



variation provides the most accurate estimation of the 
maximum slope# but with the accompanying degradation in 
estimaxes for maximum displacement and the integral of 
displacement. Again# the large errors in xhese values may 
be significantly reduced by employing a grading scheme to 
restrict the size of the larger elements as shown in Figure 
5.4 (f) . Unlike the previous problem, such a modification 

has no effect on the estimate of maximum slope because of 
the extrapolation of the element midpoint slopes, which are 
exact regardless of the element model. 
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A different version of the tapered bar problem, for 
which the displacement and strain energy criteria will nor 
produce identical refinements, involves replacing the 
concentrated tip load P with a linearly varying axially 
distributed load q(x), specified by the values at the fixed 
end and the tip . The problem may be further modified 
by reversing the bar such that the maximum slope occurs at 
the fixed end, while the maximum displacement occurs at the 




Figure 5.5 Reversed Tapered Bar with Distributed Load. 

free end as shown in Figure 5.5. The case of the linearly 
varying distributed load is included in the formulation in 
Appendix B. 

This problem was solved for a uniformly distributed load 
using the same procedure as in the previous two problems. 
The refinement models are presented in Figure 5.6 and the 
solution results in Table III. The observations are consis- 
tent with those made in the previous problems, but now one 
would likely agree that the refinement based on the strain 
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(a) Uniform 






(d) u 
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TABLE III 

Beversed Tapered Bar Solation Results 



Problem Parameters: L = 100 in 

A = 0,5 in2 

A = 10.5 in2 

E = 10x10 6 psi 

g = 100 Ib/in 



Variation 


Percentage 


Relative 


Error 


Refinement 




r*" 


Criterion 


u(max) u' 


(max) 


j,udx 


Uniform 


-5.5 


39.4 


-7. 1 


u 


-2.0 


-7.5 


-2.6 


u* 


-2.7 


-8,1 


-3.4 


U 


-3.7 


-5 . 9 


-4.7 


SED 


-11.9 


-3.4 


-15.3 


SED (mod) 


-3.1 


-3,4 


-4.0 



energy density would represent an optimal grid, provided 
that modifications are introduced to prevent any elements 
from growing excessively large. 

C. GOIDELINES FOR ONE-DISB NSIONAL GRID OPTIMIZATION 

The most important lesson to be learned from this one- 
dimensional study is that the grid optimization procedure is 
necessarily dictated by the optimization goal, or the under- 
lying purpose for performing the finite element analysis. No 
element grid can possibly provide optimum accuracy for every 
solution resultant of interest. In solving these simple 
problems, a balance has been sought for achieving adequate 
accuracy for three of the more important solution 
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the deriv- 



resaltants, with emphasis on the maximum value of -i 
ative of the dependent variable, which more often is not 
only the most important part of the solution but also the 
most difficult to obtain accurately in finite element 
analy sis. 

The important grid optimization techniques of intro- 
ducing mere degrees-cf-freedom by subdividing the elements 
or increasing their polynomial order have been intenrionally 
omitted in favor of the optimization strategy of seeking 
maximum solution accuracy for a specified number of 

degrees-of-freedom using linear elements. This is because 
such a procedure is not so straightforward in two- 

dimensional problems where the number of degrees-of-freedom 
are dependent on some geometric considerations, which do not 
appear in problems of one-dimension. Based on this choice of 
optimization strategy, it appears the strain energy density 
variation provides the most useful criterion for the adap- 
tive refinement of the initial grid. Yet all three problems 
demonstrated some pathological results that can arise when 
the elements are permitted to grow excessively large in the 
regions where the strain energy density varies the least. In 
applying a scheme to restrict the size of the largest 
elements, no mention has been made of how to determine when 
an element is excessively large. It has become the experi- 
ence of the author that any element representing over half 
of the domain should probably be considered too large, and 
measures should be employed to reszricn its size. 

It would appear, at least for these classes of problems, 
that this difficulty of decreasing accuracy of a particular 
solution parameter for successive refinements can be ignored 
by merely accepting the largest value among the cycles as 
the most accurate solution result. For example, it was 
demonstrated that the refinement based on strain energy 
density provided significant improvement in the accuracy for 




16-^ 
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ISB5 

ff 
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the maxiniuin slope but underestimated the maximum displace- 
ment even mors than the initial uniform grid, Assuming that 
the linear element model always underestimates such maxima, 
the maximum slope for refined grid and the maximum displace- 
ment for the unrefined grid could be accepted as the optimal 
results of the analysis. The fallacies of such a procedure 
are that, first, the refinement may not represent the 
optimal grid as it has been defined and, second, for self- 
adaptive finite element codes the user is provided with the 
’’optimum grid" of the final cycle and the solution results 
thereof. 

Melosh and Marcal [Ref. 21 ] have proposed an alternative 
use of the refinement criterion based on strain energy 
density variation which avoids the problem of excessive 
element growth altogether. Beginning with a reasonably 
coarse uniform grid, those elements with the greatest strain 
energy density variation are selectively refined by either 
subdividing them or increasing their polynomial order with 
the introduction of additional degrees-of-freedom. while 
such a procedure does not equi-distribute the element strain 
energy variations, it can reduce them all to some prespeci- 
fied tolerance, such as a percentage of the average element 
variations for the initial analysis. Because this procedure 
is particularly attractive for grid refinement in problems 
of higher dimensions, it will be employed extensively for 
the study of grid optimization for two-dimensional problems 
in the next chapter. 
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VI. M plication to two;:Dimehsional probLJ^s 



since investigators began working in the field of finite 
element grid optimization in the early 1970 's, nearly all of 
the effort has been devoted to the development of a system- 
atic procedure for obtaining optimal grids for two- 
dimensional problems of elasticity. Even roday rhere are 
several competing approaches to this problem and no partic- 
ular one has yet been overwhelmingly accepted as the 
preferred method of grid optimization. While it is the two- 
dimensional problem for which most of these techniques have 
been developed, their application to such can be much more 
difficult than for the one-dimensional case. Almost invari- 
ably when performing grid refinement on two-dimensional 
domains, the analyst is confronted with the problems of 
maintaining interelement compatibility and preventing severe 
element distortion. 

In selecting an appropriate two-dimensional problem for 
the application of seme grid optimization techniques and a 
comparison of their effectiveness, it is desirable that the 
test case possess the following properties: 

• the analytical solution should exist in order to 
perform reliable error analysis 

• the solution should exhibit sufficiently large 
gradients to provide a meaningful measure of the 
refinement effectiveness 

• the idealization should have one degree-of-fr eedom per 
nods and possess simple boundary conditions ~.o 
minimize the cempura rional effort involved in 
repeated solutions 
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these crireria, 



but. 



There are few problems that meet 
Saint-Venant torsion of a non-cirrular section provides a 
good test case. 

A. PROBLEM DESCRIPTION 

Consider a solid circular shaft of radius '*a” made from 
isotropic material of shear modulus G and having a circular 
groove, or keyway, of radius b along a generator of the 



Figure 6.1 Cross-section of Shaft with Keyway. 

shaft. The shaft cross-section is shown in Figure 6.1. The 
shaft is subjected to an applied torque T which produces 
an angle of twist per unit length 9. The problem may be 
solved by finding the Prandtl torsional stress function ^ 
which satisfies the governing differential equation: 




X 







(Sqn . 6.1) 
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I: 





0 cn 



subject to the Dirichlet condition that ^ = 

section boundary. The torsional stress function 
such that the shear stress at any point on the 
be expressed as: 






is defined 
domain may 






3 ^ 

y 






(Eqn. 6.2) 



For this formulation, the angle of twist 9 is prescribed, 
rather than the applied tor gue T. The torque is is calcu- 
lated from the area integral: 



T = 




(Eqn. 6.3) 



The analytical solutions of Equations 6.1 and 6.2 are 
derived by Sokolnilcoff [Ref. 22: pp. 141-143] and are 

presented in Appendix C along with the evaluation of 
Equation 6.3 and a prescribed finite element formulation. 
For this problem, the three solution resultants of interest 
for the grid optimization study are: 

(1) maximum value cf the dependent, variable, or torsion 
function \^max; 

(2) maximum value cf the gradient of the dependent vari- 
able (a quantity proportional to maximum shearing 
stress T max) ; 

(3) the area integral of the dependant variable over the 
domain (a quantity proportional to the applied 
torque T) . 

These quantities - the dependent variable, its gradient, and 
an integral thereof - are selected as representative of 
entities whose error one might wish to minimize in a finite 
element analysis. 
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B. COHPOTEB IMP LEU ENT AT ION 



As can b€ seen in Figure 6.1, the domain of this problem 
is symmetric about the x-axis, therefore the finite element 
solution need only be obtained for the upper half of the 
domain. For all of the solutions presented herein the 
problem geometry is defined by assigning the dimensionless 
ratio, b/a =3.4, and an acceptable upper limit on the anal- 
ysis cost was arbitrarily chosen to be that corresponding to 
approximately one hundred nodal points. The computation and 
assembly of the finite element matrices and solution of the 
resulting system of equations was performed using the steady 
state heat conduction operations of CAL-NPS [Ref. 23]. This 
group of subroutines comprises an efficient finite element 
code fcr solving Poisson's equation in two or three dimen- 
sions and has the additional advantage of permitting 
variable-noded isoparametric elements. 

Since there was no readily available interactive prepro- 
cessor which lent itself well to adaptive mesh refinement, 
the author had no choice but to create his own. Since the 
problem domain is simply connected, the automatic mesh 
generation was performed employing inverse mapping of a 
single cubic isoparametric element of the serendipity family 
onto the problem dcmain [Ref. 24; pp. 228-229 ]. !lappea 
boundary nodes were repositioned to conform to the actual 
domain boundary and additional nodes generated during the 
refinement process were mapped using the same procedure. 

Since the finite element code selected for this investi- 
gation provided output only for the nodal values of the 
dependent variable, it was coupled to the author's postpro- 
cessor. Such a postprocessor is necessary in the optimiza- 
tion process for computing nodal values of shear stresses 
and strain energy density, element contributions to torque 
and total strain energy, and exact results from theory. 
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c. 



ASYMPTOTIC ERROR ANALYSIS FOR ONIFORM REFINEMENT 



The concept of asymptotic convergence rate for uniformly 
refined grids was presented in Chapter 3. When the number of 
uniformly distributed deg rees-of -f reedom is sufficiently 
large, the log-log plot of the relative energy error versus 
the number of degrees-of -freedom is approximately linear in 
the asymptotic range. The slope of this line represents the 
asymptoric rate of convergence in energy. 

It so happens that relative error in rhe rorgue T of 
this problem is equal to the relative energy error and 
therefore exhibits this linear asymptotic behavior on the 
log-log plot against the number of uniformly distributed 
degrees-of- freedom. Fortunately, the other two solution 
resultants of interest behave similarly. This will prove 
very beneficial in performing the error analysis for this 
two-dimensional study for two reasons. First, because it is 
unnecessarily difficult to construct an optimal grid with 
the same number of degrees- of-freedo m as a uniform grid, the 
linear behavior of the solution resultants in the asymptotic 
range on the log-log scale permits interpolation for any 
number of degrees- of -freedo m. Then the solution results for 
a uniform grid of the identical number of degrees-of-f reedom 
provides a reference fcr comparison to determine the effec- 
tiveness of the optimization technique. Secondly, if the 
convergence cate of a particular solution resultant is 
extremely slow, as is often the case for maximum stress, it 
becomes difficult to gain an appreciation of the true effec- 
tiveness of the optimization. For example, an order of 
magnitude reduction in the relative solution error may 
require an order of magnitude increase in the number of 
degrees- of- freedom using uniform refinement, but relatively 
few additicnal degrees-of -freedom using an optimization 
technique. Therefore, it will be enlightening to extrapolate 
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(a) Linear Element Grid 




(b) Quadratic Element Grid 



Figure 6.2 Oniform Linear and Quadratic Element Grids. 
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the relative error versus d egrees-of- freedom curve to ob-air. 
a rough approximation of the number of degrees-of-f rsedom 
necessary to obtain solution accuracies similar to the 
optimal grid, but using successively finer uniform grids. Of 
course, this is only an estimation and ignores such reali- 
ties as numerical ill-conditioning and computer round-off 
error . 

A uniform grid is one for which all of rha elements are 
of the same size h and the same polynomial order p. Clearly, 
it is impossible to obtain such a grid for this particular 
domain using isoparametric mapping, but a nearly uniform 
grid may be constructed in which the elements are of approx- 
imately the same size. Such uniform grids are shown for the 
cases of linear quadrilateral elements and quadratic seren- 
dipity elements (Fig. 6.2). For this geometry, the uniform 
grid is not uniquely defined for a specified number of 
elements. This is because, in performing isoparametric 
mapping, there must be specified four keynodes on the actual 
domain boundary to correspond to the four corner nodes of 
the parenr square. Since this domain has only three 




Figure 6.3 Keynode Placement for Isoparametric Sapping. 
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vertices, the placement of the fourth keynoda is at the 
discretion of tha analyst (Fig. 6.3), and can have a notice- 
able effect on the solution results. 

The asymptotic error analysis was performed for the 
three solution resultants of interest using uniform grids of 
linear and quadratic elements. Tha results are presented in 
Figure 6.4. All of the solution resultants behaved as 
predicted with the exception of the maximum torsion function 



accuracy of this particular parameter is very strongly 
dependent upon the keyncde placement. The curve constructed 
in Figure 6.4 represents an average for several keynode 
positions. 

While Figure 6.4 is intended primarily to serve as a 
reference tool for future analyses, it provides some inter- 
esting i rforma tion: 

(1) For the cases of maximum torsion function value and 
applied torque (and energy), tha asymptotic rate of 
convergence using quadratic elements is more than 
twice that for linear elements. 

(2) While the error in torque for the quadratic case is 
always smaller than that for the linear case, the 
linear grid may provide better accuracy for the 



asymptotic range. 

(3) Both accuracy and convergence rate in the maximum 
shear stress are only minutely greater for tha quad- 
ratic element grid than for the linear one. 

However, for this last observation, the point must be made 
that for the linear element grid, the maximum shear stress 
was obtained by quadratic interpolation rather than from the 
linear shape functions. While this will greatly improve the 



value using linear elements. 




It appears that the 



maximum torsion function value 
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accuracy cf the maximum shear stress approximation, it will 
have no effect on its rate of convergence. Therefore, if 
obtaining an optimal estimate of the maximum shear was the 
purpose of the analysis, there is much to be said on behalf 
of linear elements besides their computational efficiency. 
Of course, this observation is based on uniform grid refine- 
ment, which would rarely compete favorably with the optimi- 
zation techniques to be examined. 

The reason that the rate of convergence in maximum shear 
stress is so poor using uniform refinement for this problem 




Figure 6.5 Stress Distribution on Shaft Keyway. 

can be seen in Figure 6.5. The shear stress varies greatly 
over a short distance, by increasing from zero at point A to 
its maximum value at point B. As a result, there exists a 
region of excessively large strain gradients along the 
keyway which severely hinders the rata of convergence when a 
uniform grid is employed. If the keyway radius were allowed 
to approach zero producing a singularity in the solution, it 



62 



would likely be necessary to employ even higher order 
elements via the p-version in order to achieve convergence 
using uniform refinement [Ref. 1]. 

D. PROBLEH SOLDTION WITH GRID OPTISIZATIOH 

The finite element solution of the torsion problem will 
be obtained employing the following grid optimization tech- 
niques as presented in Chapter 4: 

• Contouring 

- contours of the torsion function; linear elements 

- contours of shear stress; linear elements 

- contours of strain energy density; linear elements 

• Selective Refinement 

- h-version; linear elements 

- h-version; quadratic elements 

- p-version 

• Subdomain Isolation 

- linear elements 

- quadratic elements 

• Mesh Grading 

- linear elements 

- quadratic elements 



1 . C ontou r! nq 

The original finite element analysis was performed 
on a uniform grid of 98 linear elements, 78 nodes, and 72 
degrees-of-freedom . The finite element solution provided the 
nodal values of the torsion function ^ , from which the 
conventional nodal resultants of shear stress t , and 
applied torque T were computed. Based upon the maximum and 
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minimum values obtained for each parameter, along with 
consideration for their values along the boundary, rhe 
contours to be used for nodal placement in each case were 
selected. The number of contours for each case was chosen to 
maintain approximately the same number of degrees-of-freedom 
as for the initial analysis. 

The points for each contour value selected were 
obtained by linearly interpolating between the nodal values 
of each parameter obtained from the initial analysis. The 
contours were constructed by smoothly connecting the points 
by hand. The element layout along the contours posed the 
most formidable problem because the coarse-to- f ine tran- 
sition often resulted in severe element distortion, and it 
sometimes became necessary to degenerate quadrilateral 
elements into triangles when the transition was acute. It 
was decided that the optimal element shapes should be 
preserved along the contours in regions of highest stress. 

The contours obtained and the corresponding grid are 
presented for each of the following solution resultants; 



• torsion function 

• shear stress 

• strain energy density, SED 



(Fig. 5.6) 
(Fig. 6.7) 
(Fig. 6.8) 



The resulting grid for each of the response function 
contours produces smaller elements in the region of greatest 
stress near the bottom of the keyway and around the 
periphery of the shaft where the stress is moderately high. 
Consequently, the elements near the center where the stress 
is zero are Larger. These, of course, are the desired 
effects for an optimization criterion. A somewhat unusual 
behavior is observed at the point of intersection of the 
Iceyway and the shaft boundary where the stress is also zero. 
Apparently, the shear stress gradient is larger than the 
gradients in torsion function and the strain energy density. 
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(a) Contours 




(b) Corresponding Grid 



Contouring for tha Torsion Function. 
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Figure 6.6 




(a) Cor.tcurs 




Contouring for the Shear Stress Function, 
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Figure 6.7 



(al Conrours 




(b) Corresponding Grid 



Figure 6.8 Contouring for the SED Functi 
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resulting in smaller elements being produced in that region 
by the shear stress criterion (Fig. 6.7) than by the other 
criteria. While all of rhe grids possess to some extent the 
desirable features cf an optimal grid, the strain energy 
density function produces a far more drastic refinement 
towards the point of maximum stress, while the others repre- 
sent more moderate refinements. In fact, the SED contours 
are sc dense around the keyway that the coarse-to-f ins 
element transition scheme must include degenerate quadrilat- 
erals to avoid violating the contours. Note also thar the 
coar se- tc-f ine transition for the torsion function response 
is fairly smooth whereas the transition for the strain 
energy density and shear stress refinement tends to produce 
distorted and elongated elements. This is aggravated oy the 
fact that, unlike the torsion function, the shear stress and 
strain energy density are not monotonic over the domain. 

The solution results obtained for each grid are 
presented in the upper half of Table IV. At first glance, 
the results of the refinements are disappointing in compar- 
ison to the parenthetic values obtained using a uniform 
grid, while all three crireria produce improved accuracy for 
the maximum shear stress, the errors for the maximum torsion 
function value grow extremely large. Recalling the observa- 
tions made for xhe cne-dimensional srudy, it would follow 
that such behavior is probably due to the unusually large 
elements near the center of the domain. The entries in the 
lower half of Table IV reflecr the drastic improvement 
obtained by simply introducing a few additional degrees-of- 
freedom along those element edges which grew exceptionally 
long during the optimization process, thus increasing their 
polynomial order from one to two. Hot only did this modifi- 
cation reduce the large errors for the maximum torsion 
function estimation, but modest improvements were also 
obtained in the estimations of the other resultants as well. 
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TiBLE IV 

Contouring Solution Results 
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Once again, the selection of the optimum grid would depend 
predominantly on the optimization goal of the analysis, but 
one would likely agree that the strain energy density varia- 
tion along with some modification to restrict excessive 
element growth provides the superior refinement criterion. 

An additional word of caution is in order for the 
contouring techniques for grid optimization. Because the 
problem must be completely redefined from scratch after the 
initial analysis, the preprocessing cost can become enor- 
mous, especially if several cycles are employed to obtain 
more precise contours as some authors suggest. Unless there 
is available an interactive automatic mesh generator based 
on rhis technique, such as the one described in Reference 
[Ref. 11 ]r contouring should be abandoned in favor of some 
more easily implemented grid optimization techniques 
employing similar refinement criteria. 

2- Selective Re fine ment 

The simplest way to avoid the problems encountered 
in the contouring techniques is to perform the initial anal- 
ysis on a reasonably coarse grid and then to selectively 
refine those elements over which the strain energy density 
varies the most. The critical concern then arises as to how 
coarse the initial grid should be. If the preprocessor 
employs the necessary constraint conditions to permit the 
"father-tc-four-sons” element subdivision scheme directly, 
or if hierarchical refinement is employed, then the initial 
grid should be just coarse enough to adequately define the 
problem and to limit the number of refinement cycles neces- 
sary. The latter becomes even less of a concern if iterative 
solvers are employed. If, on the other hand, coarse-t c-f ine 
transition schemes are used to implement the h-versicn or 
only lew polynomial order elements are available in the 
p-version, then the initial grid must be sufficiently fine 
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so as net to restrict severely the amount of refinement 
which can be performed in any given cycle. Unfortunately, 
the conditions under which this investigation was conducted 
were those of the latter. 

a. The h-Version 

Selective refinement by the h-version was 
performed on both linear and quadratic element grids. For 
the linear case, the initial analysis was performed on a 
uniform grid of 55 elements, 72 nodes, and 50 degrees-of- 
freedem. The initial quadratic analysis employed an eight 
element uniform grid of 37 nodes and 20 degrees-of-freedom. 
The reason for such a grear disparity in the number of 
elements for the initial analyses is zhat subdividing a 
quadratic element introduces many more degrees-of-freedom 
than the subdivision of a linear element. These numbers were 
chosen to provide approximately the same number of dsgrees- 
of-freedem for the optimum grid of the final cycle for each 
case. The initial analysis is performed and those elements 
over which the strain energy density variation is signifi- 
cantly greater become candidates for refinement. The refine- 
ment is performed by subdividing each candidate element inro 
four new ones by constructing a coarse-to-fine transition 
zone of "buffer" elements around the refined region. 
Successive analyses and selective refinements are repeated 
until the maximum element strain energy density variation is 
approximately that of the remainder of the grid. The process 
is further improved when the nodal values of the strain 
energy density are used to indicate the general direction in 
which the refinement is to proceed. This permits multiple 
refinements in the same cycle, thereby reducing the number 
of cycles required to arrive at the optimal grid. For this 
problem, the linear grid required two refinement cycles 
while the quadratic grid required three cycles. 
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Figure 6.10 Selective Refinement Procedure - Quadratic Elements 



selective refinement process is dapicrad in Figures 6.9 and 
6.10 for the linear and quadratic element grid, 
respectively. 

The solution results for each selective rafine- 
ment cycle are presented in Table V. The most impressive 
observation to be made is the significant improvement in the 
maximum shear stress estimate for successive cycles. While 
there is also modest improvement in the accuracy of the 
torque estimate for successive cycles, when compared to the 
estimate obtained from tha uniform grid of the same number 
of degrees-of-f reedom and polynomial order, the refinement 
estimate of torque is slightly poorer. This is because addi- 
tional degrees -of- freedom are being introduced in only a 
small region of the domain but the torque, and energy, are 
global quantities. The author has no satisfactory explana- 
tion why the estimate for the maximum torsion function 
improves for successive refinements of the linear grid but 
not for the quadratic case. However, as has already been 
mentioned, this particular solution parameter appears very 
sensitive to such problem variables as nodal placemenns and 
element shapes; hence, its behavior is difficulr to predict, 
even when the refinement is applied to regions remote from 
the point where the maximum torsion function value occurs, 
as was the case in these examples. For computational 
reasons, it is desirable to restrict the number of refine- 
ment cycles to a necessary minimum. In this example, the 
quadratic grid required an additional cycle over the linear 
grid but this is because it is necessary to perform the 
initial quadratic analysis using far fewer degrees-of- 
freedom. Therefore, the early cycles of the quadratic anal- 
ysis actually represent comparatively smaller problems. 
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Selective Refinement Solution Results 
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b. The D-7ersion 



Before continuing to nhe next optimization tech- 
nique, it is worthwhile to talce a quick look at selective 
refinement employing the p-version. Because the finite 
element code used in this investigation only provided for 
element orders of one and two, the advantages of the method 
cannot be fully realized, but the effects of a single cycle 
can be examined. 

Beginning with three uniform grids with 

differing numbers of linear elements, the initial analyses 
were performed. In each case, the elements over which the 
strain energy density varied the most were transformed from 
4-noded Lagrangian elements into 8-noded serendipity 
elements by the addition of midsiie nodes. The element 

grids are shown in Figure 6.11 and the asterisks denote 
those elements for which the polynomial order was increased. 
In this example, the number of elements to be refined in 
each case was chosen so as to achieve approximately the same 
number of degrees- of -freedo m after a single cycle. 

The solution results are shown in Table VI. 
Significant improvements in the estimate of the maximum 
shear stress were achieved for each case. An improvement in 
the estimated torque was also realized for all three cases, 
but was more noticeable when the number cf refined elements 
was larger. This is because quadratic elements are far 
superior to linear elements in tae modeling of integral 
quantities, as observed in Figure 6.4. Somewhat disturbing 
is the increased error in the estimate of the maximum 
torsion function value observed in two of the three refine- 
ments even though the elements in the vicinity of its occur- 
rence were net affected. Again, this is likely attributable 
to the unusual behavior patterns of this quantity already 
discussed. 
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Figure 6.11 Selective Refinemeni: Employing p-7ersion. 
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la order to perform additional cycles of rhs 
p- 7 ersion, it would be necessary to alter the refinement 
criterion slightly. Because the element sizes do not change 
for successive cycles, the need for refinement would neces- 
sarily be based on strain energy density variation between 
nodes rather than over the elements. 

Selective refinement employing the p-version is 
most efficiently implemented hierarchically, in which case 
it acquires some attractive computational advantages. It is 
unfortunate that time did not permit further investigation 
hare, but the need for future research is evident. 

3 . Subdom a in Is olat ion 

The refinement criterion and initial procedures in 
employing subdomain isolation are identical to those used in 
selective refinement. After the candidate elements for 
refinement are identified, they are completely isolated from 
the remainder of the domain and solved as smaller subdomain 
problems. The advantages of the technique are twofold. By 
isolating the elements to be refined the solution is not 
repeated in each cycle for those elements for which the 
initial analysis is assumed adequate. Furthermore, by elimi- 
nating most of the degrees-of -f re edom over the entire 
domain, the subsequent refinement of the isolated region can 
be much greater than would otherwise be practical. 

As before, the technique was applied to both linear 
and quadratic uniform element grids. Those elements of the 
initial analysis over which the strain energy density varia- 
tion was exceptionally large were isolated to comprise the 
subdomain in each case. There were three such elements of 
the initial linear grid and two for the quadratic grid. Each 
subdomain was uniformly refined to achieve approximately the 
same number of degrees-cf- freedom as the initial analyses. 
The process is depicted in Figure S.12 for the linear grid 
and Figure 6.13 for the quadratic case. 
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(b) Rsfinsmant of Isolarsd Subdoma 



Figure 6.12 Subdomain Isolation - Linear 
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(a) Initial 3rid 




Sabdonain Isolation Quadratic 



8 1 



Figure 6.13 
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In performing subdomain 


isolation 


the 


gcve rninq 


equations remain the same, while 


only the 


doma in 


and the 


boundary conditions are altered. 


'When 


the 


subdomain 


boundary has nodes common to the 


initial g 


rid. 


then the 



boundary values for those nodes are simply the solution 
values obtained from the initial analysis. The boundary 
values arising from the introduction of new boundary nodes 
during rhe refinement process musr be generated by interpo- 
lation of the solution results of the initial analysis. One 
of the options for an interpolation scheme is simply to use 
the shape functions of the unrefined elements. This may not 
be desirable in the case of linear elements, sc a higher 
order interpolation may be employed. In this example a third 
order Lagrangian polynomial was convenient for the linear 
case since there are four nodes from the initial analysis 
along the right-hand boundary of the subdomain. Since there 
are only two such nodes on the upper boundary, it is neces- 
sary to "borrow" some adjacent nodal values from the 
discarded portion of the domain in the generation of new 
boundary values using higher order interpolation. 

The solution results presented in Table VII are 

quite remarkable. In a single cycle, the solution accuracy 
for the maximum shear stress has increased by a full order 
of magnitude. No other optimization technique examined in 
this investigation produced such improvements. Note that 
the higher order polynomial interpolation for the boundary 
values did improve the solution results for the linear case. 
One of the disadvantages of this technique is that the 

refinement can produce no improvement in the estimation of 
local quantities outside the subdomain. As in this example, 
the estimation of the maximum torsion function value 

obtained from the initial analysis must be accepted as the 

optimum since it occurs outside the subdomain. Furthermore, 
since its value is predominantly affected by refinements in 
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regions of high gradients, ir is doubtful that isolating a 
new subdomain using the elements adjacent to its point of 
occurrence would noticeably improve its accuracy. However, 
since the torgue is a globally computed quantity, refinement 
will improve the accuracy of its contribution from the 
subdomain resulting in improvement of its accuracy overall 
as observed in Table VII. It is this strictly local nature 
of the subdomain isolation technique which restricts its 
applicability, 3ut if the optimization goals are well 
defined and it is understood under which conditions and for 
what parameters it is effective, it can be an extremely 
powerful grid optimization technique. 

^ • Ma sh 3 r a di n g 

while mesh grading is nearly always performed on an 
”a-priori" basis, it may also be employed adaptively to 
provide a simple grid optimization technique. After an 
initial solution has been obtained, the analysis may be 
repeated using various combinations of grading ratios in 
order to achieve a more uniform distribution of the element 
strain energy density variations. Here the grading ratio 
refers to the constant ratio of adjacent element lengths 
along a boundary of the domain to which grading is applied. 
There are several drawbacks to the technique, the first 
being that a good combination of grading ratios may be 
difficult to obtain in a reasonable number of cycles. The 
other difficulty is that if smaller elements are desired in 
more than one region the domain must be refined and 
constructed by subregions. 

Unfortunately, the domain of this problem is not 
well suited for mesh grading since it possesses no favorable 
geometric symmetry. Hence, the resulting element elongation 
and distortion would become severe for larger grading 
ratios. For simplicity, the nodal placements will be biased 



84 



only along ths two domain boundarias adjacent to the point 
of maximum stress and the same grading ratio will be applied 
for each. This will result in small elements near the bottom 
of the keyway and large elements along the periphery of the 
shaft. While this is not the most desirable grid topology, 
it will produce a mors uniform distribution of the element 
SED variations. 

The technique was applied to both linear and quad- 
ratic element grids starting with a uniform mesh and succes- 
sively increasing the grading ratio until the elements along 
the shaft periphery exhibited SED variations as large as 
those for the elements along the key way. In both cases this 
condition occurred beyond the point where excessive element 
elongation would be expected to produce numerical inaccura- 
cies. Graded meshes for selected grading ratios r, are shown 
in Figure 6.14 for linear elements and Figure 6.15 for quad- 
ratic elements. The solution results are presented in 
Table VIII. As can be observed, the maximum shear stress 
estimate improves for each successive increase in the 
grading ratio. However, the cost of such improvement is the 
accompanying degradation in the estimate of the maximum 
torsion function value. This is to be expected since the two 
maxima occur at different locations in the domain and there- 
fore decreasing the size of the elements in the vicinity of 
one will necessarily increase the element size near the 
other. Note that this degradation is not nearly so severe 
for the quadratic elements, and the accuracy actually 
improves for a low value of the grading ratio. This is 
because the higher order interpolation can better accomodate 
larger elements. For both linear and quadratic element grids 
the optimal accuracy in torque estimation occurs for the 
moderately graded meshes. 
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Figure 6.14 Mesh Grading - Linear Elements 
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Figure 6.15 Mesh Grading - Quadratic Elements 



TABLE VIII 

Mesh Grading Solution Results 





Percentage Relative Error 


Grading 

Ratio 


^max Tmax /G9 I / G0 



Linear Element Grids ( 78 elements, 98 nodes, 

72 degrees-of-freedom ) 



1.0 


0.060 


-6.06 


-1. 77 


1. 1 


0. 16 1 


-2.63 


-1 . 56 


1.2 


0.389 


-1.03 


-1 . 77 


1.3 


0.679 


-0.484 


-2. 20 



Quadratic 


Element Grids 


( 28 elements, 
76 degrees-of 


107 nodes 
-freedom 


1.0 


-0.0093 


-5. 26 


-0.0116 


1 . 1 


0.0063 


-1.85 


-0.0064 


1.2 


0.0162 


-0.606 


-0. 0091 


1.3 


-0.0246 


-0.413 


-0.0179 



It is likely that some of the error is attributable 
to the numerical inaccuracies due to element distortion. 
When applying a grading technique the analyst should seek an 
equitable balance between the refinement criterion and the 
grid topology. 
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Finally, since the grading ratio is usually applied 
to the nodal separation rather than the element sdge 
lengths, it is advisable to reposition the midside nodes of 
quadratic elements so that they lie near the center of the 
element edges. This will generally improve the accuracy of 
all the solution parameters, especially if the grading is 
somewhat extreme. 

B. GOIDELINES FOB T aO-DIME NSIONAL 3BID OPTIM IZiTIOH 

In order to provide some guidelines for obtaining 
optimal finite element solutions for two-dimensional prob- 
lems it is helpful to compare the solution results obtained 
for this problem employing the optimization techniques 
available. Such a comparison is presented in Table IX. The 
upper portion is for those techniques for which t-he inirial 
analysis was performed using linear elements and the lower 
portion using quadratic elements. Note that all of the grids 
employ approximately the same number of degrees-of-freedom, 
which was the chosen measure of analysis cost in this 
investigation. 

In making this comparison it is important to understand 
just how significant a change in error actually is. If the 
convergence rate of a solution parameter is very slow, even 
a small reduction in the error may require many mors 
degrees-of-freedom. For this reason, the numbers in paren- 
theses have been included by each error to provide a rough 
approximation of the number of degrees-of-freedom required 
to obtain similar accuracy using a uniform grid of the sane 
number of degrees-of-freedom and elements of the same poly- 
nomial order as the initial analysis. In some cases, the 
analyses were not actually performed but instead the numbers 
in parentheses were obtained by extrapolation of the error 
versus degrees-of-freedom curve (Fig. 6.4), and ignoring 
round-off and ill-conditioning effecrs. 
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The first observation to be maie from Table IX is that 
while all of the optimization techniques produced signifi- 
cant improvemenxs in the accuracy of the maximum derivative 



might even conclude that grid optimization is not cost 
effective in the computation of these values since the 
uniform grid provides estimates which are nearly as accu- 
rate, and in some cases better, than the optimum grid. This 
conclusion would be correct if the solution maximum and its 
integral were the only resultants of interest in performing 
the analysis. Since the purpose of this study was to find an 
optimum grid which produced acceptable errors for all the 
resultants, the uniform grid is clearly inadequate. 
Moreover, since in the majority of engineering problems it 
is the derivative of the solution variable which is of 
primary interest, it deserves special consideration in 
making this comparison. 

Furthermore, one might conclude that the reason the 
error in the aajcimum solution variable is larger for the 
optimal grid is because the strain energy density variation 
criterion always concentrates the degrees-of -freedom in the 
vicinity of the maximum derivative value, which in this case 
does not coincide with that of the maximum solution variable 
value. However, such a conclusion is incorrect and might 
erroneously lead one to attempt refinement where the maximum 
solution variable occurs in an effort ro improve its esti- 
mate. Other investigations have revealed that in nearly all 
cases the maximum accuracies for all three solution resul- 
tants are obtained by refinement in the regions of highest 
strain energy density variation [Bef. 11]- The more likely 
source of increasing error for the optimal grids is rhe 
element distortion which was encountered in all bur two 
techniques, selective p-version refinement and subdomain 



quantity Tmax, 
solution quantity 



the same cannot be said for the maximum 
and the integral quantity T. One 
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isolation. Sach distortion can be avoided but would require 
more sophisticated refinement techniques than were available 
for this investigation. 

A reasonable choice for the optimum grid in Table IX 
would be one for which all three values in parentheses are 
as large as reasonably possible talcing into consideration 
the number of cycles required to provide such accuracy. 
Based on such a criterion, the author is partial to subdo- 
main isolation for the solution of two-dimensional problems 
using linear elements, and selective refinement for finite 
element solutions using quadratic elements. Clearly, before 
a concrete recommendation could be made for a wide range of 
applications, many more problems would have to be studied, 
but these two techniques were fairly simple to implement for 
a standard finite element code and the accelerated conver- 
gence of the solution resultants of interest was impressive. 
Conceivably, even greater solution accuracies might be 

obtained by using two or more of these techniques in 
combination . 

Here again the crucial element in selecting the proper 



optimization 


strategy is the 


precise 


definition of 


the 


purpose for 


which the finite 


element 


analysis is to 


be 


performed. 


The results of Table IX t 


end to support 


the 



following recommendations and conclusions; 

(1) Regardless of the optimization strategy chosen, 
higher order elements are indispensable if high 
accuracies for integral solution quantities are 
desired. 

(2) If the maximum derivative of the solution variable 
is of greatest concern, the strictly local refine- 
ments employing subdomain isolation techniques can 
provide exceptional accuracy for a minimum number of 
refinement cycles. 
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(3) If the Baxiaum solution variable value occurs ar a 
point in the domain removed from the vicinity of the 
maximum derivative value, then its best estimate 
will likely be obtained using a reasonably fine 
uniform grid and selectively subdividing elements in 
the regions cf large strain energy density 
variations. 
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VII. co Nciosioa MS recqmmb ndat ions 



The purpose of this paper has been to present an over- 
view of some readily employed finite alemenr grid optimiza- 
tion merhods and zo demonstrate their effectiveness in the 
application to some simple problems. This work is by no 
means ail inclusive and the subject is still in its infancy. 
While there are many competing approaches to the problem, 
there is much more research to be done before any one 
becomes widely accepted as a standard analysis tool. Because 
of the limited time and resources available, some of the 
more sophisticated refinement criteria and techniques which 
have been developed have not been examined in any detail. 
Instead, the approach has been to examine those techniques 
which can be easily incorporated in a basic finite element 
code. However, it is likely that some recently developed and 
rather elaborate self-adaptive grid optimization codes will 
soon be available. 

Also, this paper has not examined the important classes 
of problems in dynamic and nonlinear analysis. There is 
considerable ongoing research in the extension of these 
techniques to such problems, but the increased complexity is 
evident. For example, in vibration analysis there is an 
optimum grid for each unique eigenvalue, but it is for these 
types of problems that grid optimization is most promising. 

At the beginning of this paper it was stated that the 
goc,l of grid optimization was to obtain maximum solution 
accuracy for a given analysis cost. Throughout this paper it 
has been shown that, prior to successfully embarking upon 
such a strategy, the underlying purpose of the analysis must 
be explicitly defined. Hopefully, it has been demonstrated 
that grid optimization is by no means an unrealistic goal 
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and is far mDre attractive than the non-adaptive practices 
widely used today. 

The following are recommendations for future research 
topics: 

(1) Investigation of more sophisticated refinemenr 

criteria based on element residuals and reliable 
error estimates. 

(2) Investigation of grid optimization rechnigues 

employing adaptive application of the p-version, 

(3) Implementation of a finite element preprocessor for 
performing hierarchical grid refinement. 

(4) Implementation of a self-adapzive finite element 
code. 

(5) Application of grid optimization techniques to prob- 
lems of dynamic and nonlinear analysis. 
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4PPEHDIX A 

FORMOLATION OF THE ELASTIC CABLE PBOBLEM 
A. PBOBLEM STATEMENT 

Consider a perfectly elastic cable initially stretched 
between two fixed points a distance 2L apart and under 
tension T. If the cable bears a distributed load per unit 




Figure A. 1 Tension Cable Under Distributed Loading. 



length f (x) as shown in Figure A.1, the governing differen- 
tial equation for the downward defection v(x) is: 



T 




+ f (X) = 0 



(Egn. A.1) 



subject to the essential boundary condition: 



V (X = ±L) = 0 



(Egn. A. 2) 
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If th€ distributed load is a supportive load provided by 
a Hinkler foundation of modulus k such that 

f(x) = - k v(x) 

and if a concentrated load P is applied at the midspan. 
Equation A. 1 becomes: 




(Eqn. A. 3) 



subjecr to the natural boundary condition: 




(Eqn. A.U) 



B. PROBLEM SOLOTION 

The analytical solution of the two-point boundary value 
prcblem is: 



V (X) 



- P/2 
" TA 



[tanhAL coshAx - sinhAx] 



where 



= k/T 



(0 < X < L) 
(Eqn. A. 5) 



The finite element solution is obtained by the Galerkin 
formulation using the consistent rarher than the lumped 
approximation for the distributed loading. 
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appendix b 

FOEMOLATION OF THE TAPERED BAR PROBLEM 
A. PROBLEM STATEMENT 

Consider a tapered bar of length L and consnant modulus 

of elasricity E fixed at one end. The cross-sect ional area 

A(x) varies linearly from A at the fixed end to A. at the 

tip. Let the bar be loaded axially by a concentrated tip 

load P, and a distributed load for which the intensity g (x) 

varies linearly from g at the fixed end to g at the tio as 

o t ' 




Figure B.1 Tapered Bar with Applied Loads. 



shewn in Figure B. 1. The governing differential eguation 
for the axial displacement u(x) is; 



2 i 



A ^ 
^ ox 



+ a = 0 



(0 < X < L) 



(Sqn. 3.1) 
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subject to the essential boundary condition: 



u (X = 0) = 0 



(Sgn. B.2) 



and the natural boundary condition: 



^ P 

^ x=L " EAt 



(Eqn. B.3) 



B. PEOBLBH SOLUTION 

Let a = 1 - A /A 

t o 



a nd 



6 = 1 - gt/q 



O* 



ri 

For F (X) = P ♦ I q(x) dx the solution is: 

J X 



u(x) =- 



X 
0 



F(x ) 

A(x) 



dx 



aA E 
o 



^{1- - L(l- - 

2 a q 



z„a- 2f) * ( 1 - 

(0 <_ X <_ L)'' 
(Eqn. B.4) 



S 



6x^ 






The finite element solution is obtained by the Galerkin 
formulation using the consistent rather than the lumped 
approximation for the distributed loading. 
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AP PENDIX C 

PORMOLATION OP THE TORSION PROBLEM 
A. PROBLEM STATEMENT 

Consider a solid circular shaft of radius "a" and shear 
modulus G, with a circular groove, or keyvay, of radius b 
along a generator of the shaft with the cross-section shown 



y 




Figure C. 1 Cross-section of Shaft with Keyway. 

in Figure C. 1. An applied torque I will produce an angle of 
twist per unit length 9. The equilibrium condition is 
satisfied if a torsional stress function i{; exists such that 
the shear stress components are: 



X 

ZX 




and 



T 

zy 



-G9 



3 X 
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The governing differential equation for the torsional stress 
function rp (x,y) is: 



subject to the Dirichlet condition, ^ = 0 on the boundary. 
B. PHOBLEM SDLOTION 

The solution of Equation C.1 [Baf. 22; pp, 141-143] is; 

2 

P(x,y) = a(x-b^ — - — ) - h(x^+y“) + ^ (Sqn. C.2) 

x2+y2 2 

The maximum shear stress oc urs at point A (Fig. C.1) and is: 



9^ 



0 



(Eqn. C.1) 






Tnax = G0 (2a - b) 



(Eqn. C . 3 ) 



The applied torque computed from the area integral: 




(4a‘*-8a2b2-2b'*)a + (2a^h+7ah^ )-Un a 



(Eqn C.4) 



where a = arcos (b/2a) . 

The strain energy per unit length of shaft is: 



U» 




(Eqn. C.5) 
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The variational formulation of the 
approximation is presented in detail in 
Reference 25. 



finite element 
Chapter 6 of 
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